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Abstract: We study the general relativistic non- linear dynamics of self-gravitating ir- 
rotational dust in a cosmological setting, adopting the comoving and synchronous gauge, 
where all the equations can be written in terms of the metric tensor of spatial hyper-surfaces 
orthogonal to the fluid flow. Performing an expansion in inverse powers of the speed of 
light, we obtain the post-Newtonian equations, which yield the lowest-order relativistic 
effects arising during the non-linear evolution. We then specialize our analysis to globally 
plane-parallel configurations, i.e. to the case where the initial perturbation field depends 
on a single coordinate. The leading order of our expansion, corresponding to the "Newto- 
nian background" , is the Zel'dovich approximation, which, for plane-parallel perturbations 
in the Newtonian limit, represents an exact solution. This allows us to find the exact 
analytical form for the post-Newtonian metric, thereby providing the post-Newtonian ex- 
tension of the Zel'dovich solution: this accounts for some relativistic effects, such as the 
non-Gaussianity of primordial perturbations. An application of our solution in the context 
of the back-reaction proposal is eventually given, providing a post-Newtonian estimation 
of kinematical back-reaction, mean spatial curvature and average scale-factor. 
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1. Introduction 

Large-scale structures in the Universe grew by gravitational instability around primordial 
seed perturbations generated during inflation. In this context, it is assumed that the 
matter mainly consists of Cold Dark Matter (CDM), which, prior to caustic formation, 
behaves like a perfect fluid of dust. The study of the gravitational dynamics of CDM can 
be performed by different techniques and approximations, depending on the specific range 
of scales under investigation. 

The most widely used approximation is standard perturbation theory (see e.g. Ref. Q), 
which is an expansion in powers of the amplitude of fluctuations around a Friedmann- 
Robertson- Walker (FRW) background. The range of applicability of this perturbation 
technique is that of small fluctuations around a FRW background, with no extra limitations 
on scale. Going to higher perturbative orders generally helps to follow the gravitational 
instability process on a longer time-scale and to account for non-linear phenomena, i.e. to 
resolve smaller scales. 

Another important approximation scheme is the Newtonian one, which is known to 
produce accurate results as long the gravitational field is weak and the motion of both 
sources (i.e. CDM fluid elements) and test particles is slow. In the cosmological framework. 
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the range of validity of the Newtonian approximation concerns scales much larger than the 
Schwarzschild radius of collapsed bodies and much smaller than the Hubble horizon (see 
e.g. Ref. 0). This approximation is usually carried out in the Poisson gauge, where 
the FRW metric is perturbed by a single lapse function which represents the peculiar 
gravitational potential sourced by matter density fluctuations. This approach provides 
the well-known Newtonian Eulerian description of cosmological structure formation in the 
expanding universe. The extension of this approximation to deal with fast moving test 
particles (like photons) is the so-called weak-field limit, where a scalar potential is accounted 
for in the spatial part of the line-element. 

In the context of Newtonian gravity, departures from the FRW background may be also 
analyzed in the Lagrangian frame. In the unperturbed situation, observers comoving with 
the Hubble flow are Lagrangian, but when density fluctuations are taken into account, they 
become Eulerian observers with respect to the peculiar motion. The equations describing 
the dynamics of CDM, which we model here through irrotational dust, are written in La- 
grangian coordinates, in terms of the displacement vector S from Lagrangian to Eulerian 
coordinates. As in the Eulerian case, it is impossible to work out the general analytical 
solution for S and a perturbative approach is again introduced by means of an expansion 
in powers of the displacement vector, the background being once more represented by the 
FRW model. The linear result is the so-called Zel'dovich approximation [^. The peculiar- 
ity of this treatment, at any order, is that, while the displacement vector is calculated from 
the equations at the required perturbative order, all the other dynamical variables, such as 
the mass density, are calculated exactly from their non-perturbative definition. Since the 
equations in Lagrangian coordinates are intrinsically non-linear in the density, what comes 
out is a fully non-linear description of the system, which, though not being generally exact, 
"mimics" the true non-linear behavior. This perturbation treatment basically exploits the 
advantages of the Lagrangian picture, leading, in particular, to a more accurate descrip- 
tion of high density regions. Its limitations are generally set by the emerging of caustic 
singularities. 

The Newtonian analysis can be improved by a post-Newtonian (FN) approach. In the 
Newtonian approximation two conditions must be satisfied: weak gravitational field and 
slow motion. Relaxing one of these two conditions leads to two different improvements of 
the Newtonian approximation. The FN approximation is suitable for a system of slowly 
moving particles bound together by gravitational forces. This situation does not imply weak 
gravitational fields. Therefore, in cosmology, the FN approximation can be used to account 
for the moderately non-linear field generated during the non-linear stage of the evolution of 
matter fluctuations. The FN approximation is formally obtained via an expansion in inverse 
powers of the speed of light c. The FN approximation in the cosmological framework has 
been studied in Eulerian (or more generally non-comoving) coordinates in Refs. [Q, ^, 
while in Ref. Q a hybrid approximation scheme is proposed which upgrades the weak-field 
limit of Einstein's field equations to account for post-Newtonian scalar and vector metric 
perturbations and for leading-order source terms of gravitational waves, while including also 
the first and second-order perturbative approximations. An alternative approach has been 
recently proposed in Ref. |To[| , whose authors introduce an effective fluid description for 
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small-scale non-linearities in the framework of General Relativity. The PN approximation 
to cosmological structure formation within the Lagrangian picture has been studied by 
Matarrese &: Terranova in Ref. |11], whose approach will be the basis of our analysis here. 

The aim of this paper is to find a PN metric describing the non-linear stages of the 
gravitational instability in an Einstein-de Sitter universe, adopting the Lagrangian picture 
i.e. the synchronous and comoving gauge. We choose the leading order solution of our 
expansion to be the Newtonian metric related to the Zel'dovich approximation, in the very 
special case of globally planar dynamics due to a perturbation seed depending on a single 
Lagrangian coordinate. 



Our approach here is very different from those of Refs. [|T^, 13|, who proposed rela- 
tivistic generalizations of the Zel'dovich approximation. Ref. [^] introduced a relativis- 
tic, tetrad-based, perturbative approach, which is then solved to linear order and used 
to obtain non-perturbative expressions for the velocity-gradient tensor and mass density. 
The solution of Ref. ||l^ is instead equivalent to a relativistic second-order perturbation 
theory treatment in the synchronous and comoving gauge, in which all quantities (metric, 
velocity-gradient tensor and mass density) are calculated at second order, thereby partially 
missing the non-perturbative character of the Zel'dovich approximation. Our approach in- 
stead aims at obtaining a non-perturbative description of both metric and fluid properties 
(velocity-gradient tensor and mass density), within the post-Newtonian approximation of 
General Relativity. In our approximation scheme the Zel'dovich solution represents the 
Newtonian background over which PN corrections can be computed as small perturbations. 

The plan of the paper is as follows. In Section 2 we introduce the relativistic equations 
governing the dynamics of irrotational dust. Section 3 deals with the Newtonian limit. 
Section 4 is devoted to the PN approximation: we first describe the Zel'dovich Newtonian 
background chosen as a basis for our PN expansion, showing explicitly all the dynamical 
quantities involved. Then, we obtain and solve analytically the PN equations for the metric. 
In particular, we show that our solution accounts for the non-Gaussianity of primordial 
cosmological perturbations. To conclude this section, we briefly discuss the accuracy of 
our perturbative expansion. Section 5 deals with an application of our solution: we use 
our metric in the context of the back-reaction proposal. To be more specific, we give an 
estimate of the PN contribution to the kinematical back-reaction, mean spatial curvature 
and average scale-factor, as defined in Ref. W^- Section 6 provides a comparison between 



our PN solution and the Szekeres metric. Conclusions are drawn in Section 7. 



2. Relativistic dynamics of irrotational dust in the Lagrangian picture 

In this section we recall the relativistic equations governing the evolution of irrotational 
dust following the Arnowitt-Deser-Misner (ADM) formalism [^] . Calculations are made in 
the synchronous and comoving gauge. In this gauge both the world-lines of the observers 
coincide with the geodesies of the matter and the proper time along the geodesies coincides 
with the cosmic time. This reference frame is therefore Lagrangian. Actually, the possibility 
of making these two gauge choices simultaneously is a peculiarity of irrotational dust, which 
holds also in the non-linear stage. 
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The line element is given by 

ds^ = -c^dt^ + ha0{t, q)d9" V , (2-1) 

where qa are the spatial Lagrangian coordinates and t is the proper time of our fluid 
elements. The matter stress-energy tensor reads T'^'^ = q(?u'^u^ with q the mass density 
and n° the fluid four-velocity (normalized to u^Ua = —1), with components = (1, 0, 0, 0). 
The ADM equations can be written in terms of the velocity-gradient tensor defined as 

ei = cu".f, = ^h^^'K^ , (2.2) 

where dots denote partial diff'erentiation with respect to the time t and semi-colons stand 
for covariant differentiation. The tensor 0^ represents the extrinsic curvature of constant 
t hyper-surfaces. 

The 10 Einstein equations split into 4 constraints and 6 evolution equations. The 
time-time component of the Einstein equations is the so-called energy constraint of the 
ADM approach, which reads 

- Q'^Q^ + (% = WttGq , (2.3) 

where the volume-expansion scalar B is the trace of the velocity-gradient tensor and ^^^R is 
the trace of the 3-dimensional Ricci tensor ^^^R'^ of the spatial hyper-surfaces of constant 
time. 

The space-time components give the momentum constraint, 

Q%a = 0,/3 • (2.4) 

The space-space components govern the evolution of the velocity-gradient tensor and rep- 
resent the only truly evolution equations. They read 

+ 0G^ + = AttGqS'^ . (2.5) 

Taking the trace of the latter equation and combining it with the energy constraint yields 
the Raychaudhuri equation for the volume- expansion scalar 0, 

+ 2S2 + l02 + 47rG^ = O, (2.6) 
o 

where = G^- (1/3)65^ is the s/iear tensor and = (1/2)E^S^ its magnitude. Finally, 
the continuity equation for the stress-energy tensor, T^p = gUaUp, leads to 

g = -eg, (2.7) 

whose solution reads 

h being the determinant of the spatial metric. Here and in what follows quantities with a 
subscript "in" are meant to be evaluated at some initial time tin- 
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2.1 Conformal rescaling and FRW background subtraction 

With tiie purpose of studying gravitational instability in the Einstein-de Sitter background, 
we first factor out the homogeneous and isotropic expansion of the Universe. To this aim, 
we perform a conformal rescaling of the metric with the conformal factor a(T) oc r^, the 
scale factor of the Einstein de-Sitter model. This procedure involves the transformation of 
the proper time t to conformal time r via 

dr = ^ (2.9) 
a{t) 

and it is discussed in greater detail in Ref. The line-element is written in the form 

ds"^ = a2(r) \-c^dT^ + 7„^(r, q)dq"dq'^] , (2.10) 

where 7q-^(t, q) = {t (r) , q) /a? {t (r)) is the conformal spatial metric. 

In order to factor out the Einstein-de Sitter background, the isotropic Hubble flow, with 
velocity-gradient tensor (a' /a)6'^, is subtracted from the velocity-gradient tensor, leading 
to a tensor which describes gradients of the peculiar velocity only: 

= - = ^T^^^T^/? , (2.11) 

where u°' = (1/a, 0,0, 0) and primes denote partial differentiation with respect to the 
conformal time r. The tensor represents the extrinsic curvature of constant r spatial 
hypersurfaces. The matter content is also written in terms of the density contrast, defined 
as the dimensionless deviation of the matter density from that of the Einstein-de Sitter 
background, £>b(r) = 3/(27rGa^T^), 

^(r,q) - (2.12) 

0b{r) 

The ADM equations are then recast in a more convenient form describing the evolution of 
the peculiar velocity-gradient tensor. The energy and momentum constraints become 

Q 24 
- + + c27^ = — 5 (2.13) 

= , (2.14) 

where Va is the covariant derivative associated with 

is the conformal three-dimensional (scalar) curvature. The evolution equation becomes 

+ + m'^ + + c^TZ'^ = ^66'^ . (2.15) 
After replacing the density from the energy constraint, this equation is written as 

= 7?°' + -7?^ + ^^J^ + 1 (iJM^j; _ S<^ + 



T P P 4 

+^[47^^-7^5^]. (2.16) 



c2 
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The Raychaudhuri equation describes the evolution of the pecuhar volume expansion scalar 
and reads 

^' + _^ + ^/;^^ + _j = 0. (2.17) 



Finally, the solution of the continuity equation ( |2.8| ) is written in terms of the contrast 6 
as 

5(q,t) = (l + Mq))^^^-l, (2.18) 

where 7 is the determinant of the conformal spatial metric 7^/3. 

The main advantage of this formalism is that there is only one dimensionless variable in 
the equations, namely the spatial metric tensor 7^^, which is present with its partial time 
derivatives through and with its spatial gradients through the spatial Ricci tensor 7^^. 
A relevant advantage of having a single tensorial variable, for the following PN expansion, 
is that there can be no extra powers of c hidden in the definition of different quantities. 

3. Newtonian approximation 

The Newtonian equations in Lagrangian form can be obtained from the full relativistic 
equations of the previous section by an expansion in inverse powers of the speed of light. 
As a consequence of our gauge choice, however, no odd powers of c appear in the equations, 
so the expansion parameter is 1/c^. The spatial metric is then expanded in the form 

laP=lal3 + -^Wal3 + 0(^^^ . (3.1) 

Let us first concentrate on the Newtonian metric 7q,/3. The Newtonian limit (c — )• oo) 
of the energy constraint and evolution equation requires that the spatial Ricci tensor is 
zero. Thus, the Newtonian result of vanishing spatial curvature is recovered (|jl^, [p^]). 
This important conclusion implies that 7^^ can be transformed to the Euclidean metric 
6ai3 globally. In other words, at each time r there exist global Eulerian observers comoving 
with the Hubble flow for which the components of the metric are Sa^. This means that, 
according to the tensor transformation law, we can write 

7q/3 = Sf^uJ^Jp , (3.2) 

where 

is the Jacobian of the transformation 

x(q,r) =q + 5(q,T) (3.4) 

and x^^, (/i = 1,2,3), are the Eulerian (comoving) coordinates. The tensor dS^/dq°' is 
called the deformation tensor. All the information about the motion is contained in the 
mapping between the comoving Eulerian and Lagrangian coordinates, whose evolution is 
described by the Raychaudhuri equation and the momentum constraint. In fact, contrary 
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to the evolution equation and the energy constraint, the Raychaudhuri equation and the 
momentum constraint contain no exphcit powers of c, thus preserve their form in going to 
the Newtonian hmit. These equations determine the background Newtonian metric 7^^: 
they are indeed the Newtonian equations in the Lagrangian approach to irrotational dust 



dynamics |16]. Thus, it makes sense to recast these equations in terms of the map between 
comoving Eulerian and Lagrangian coordinates. One can equivalently write the equations 
in terms of either the displacement vector |17] and deformation tensor or in terms of the 



Jacobian matrix of the map |11]. Here, following Ref. [11|, we consider the Newtonian 



equations in terms of the Jacobian matrix ( p. 3D . They read 

2 ^ _6_ / 1 
T J tA J 



J^J!i = -,{!--] (3.5) 



e"^^JjJ-;^ = 0, (3.6) 

where is the determinant of the Jacobian matrix. We assumed for simplicity 5in = and 
used the residual gauge freedom of the synchronous and comoving gauge to set = 1 in 
the Newtonian limit |11|. 

It is impossible to follow the dynamics in an exact analytical way and one is forced to 
use approximation techniques to solve these equations. Here the perturbative expansion is 
performed in terms of the Jacobian matrix. The key point is that a slight perturbation of 
the Lagrangian particle paths carries information on the non-linear dynamics, as may be 
seen directly from the expression of the density contrast 

5 = -^-!, (3.7) 

which that, even if the displacement vector is small, the corresponding density contrast 
could be large. In other words, already the first-order Lagrangian equations describe the 
mildly non-linear regime of the gravitational instability. 

Expanding and solving the previous equations to first order, the well-known Zel'dovich 
approximation is recovered: 

x(q,r) = q + I?(r)V$(q), (3.8) 

where D{t) = t'^/t^^ is the growing mode solution for the Einstein-de Sitter model and 
<I>(q) has to be ascribed to the chosen initial conditions. 



4. Post-Newtonian approximation 

4.1 Characterization of the Newtonian background 

The starting point of our PN expansion is the Newtonian background described by the 
Zel'dovich approximation, with the peculiar gravitational potential depending on the con- 
formal time and on the Lagrangian coordinate qi only. As it is well known, in the particular 
case of planar perturbations the Zel'dovich approximation yields an exact solution of the 
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Newtonian equations. The Zel'dovich solution is given by 

xi = qi+ D{T)di'^in{qi) 
X2 = q2 

X3 = qs- (4.1) 

where di = d/dqi and the potential <I>in is defined so that df^in = —^in and is related to 
the initial peculiar gravitational potential <p by the cosmological Poisson equation, yielding 



The Jacobian matrix is 



$in = --^7^ . (4.2) 

4:TTGaf^Qb;m 

/ 1 - T^dfif/G \ 
J^=i 1 j ■ (4.3) 

For the following calculations, it is useful to define the function / 

/ := D{T)dl^,^ = -^^iV (4.4) 

and the function rj = ln(l + /). Hereafter, the peculiar gravitational potential is meant 
to be evaluated at the initial time Tin and the subscript "in" is dropped for notational 
convenience. 

The components of the metric in Lagrangian coordinates are given by (recall Eq. ( |3.2| )) 

7„/3 = K + D{T)d'^d^<!>] [<5^ + D{T)d^dp^] , (4.5) 

where we used the fact that at first order in the displacement vector covariant and partial 
derivatives with respect to the coordinates qa coincide, since the Christoffel symbols are 
second-order quantities. In our case $in = ^miqi) and for the Zel'dovich metric we find 

/(i-T2aiV/6)2 \ 

7a/3 = I 1 I (4.6) 

or in more compact form 

/(! + /)' \ 

7a/3 = I 1 1 • (4.7) 

It is important to keep in mind that this Newtonian metric is non-linear with respect 
to the peculiar gravitational potential, thus it characterizes the mildly non-linear stage of 
the gravitational instability. Starting from this metric at first order in the displacement 
vector, all the other dynamical variables are calculated exactly. The only non-vanishing 
component of the peculiar velocity-gradient tensor is 'di = r( and for the shear tensor we 
have a\ = 27/73 and 'o\='a\ = —a\/2. 
Finally the density contrast (|3.7| ) takes the form 

7=Y^-l. (4.8) 
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4.2 Post-Newtonian expansion 

For the PN expansion, we write the metric in the form 

1 

lap = lal3 + ■ 



(4.9) 



For consistency with our Newtonian background solution, in which the pecuhar gravita- 
tional potential depends only on qi, the PN perturbation Wap can be assumed to depend 
on the conformal time and on the Lagrangian coordinate qi only. 

A clue for the form of the perturbation w^p follows from the initial conditions of 
cosmological perturbations. Even though our case of globally planar dynamics is purely a 
toy-model we prefer to assume that our initial perturbation seed is consistent with having 
been generated during inflation in the early Universe, so that our analytical results can 
give us a hint of what happens in the real Universe, where the initial perturbation seed 
is a random field that depends on all the spatial coordinates. We then set our initial 
conditions at the end of inflation, effectively coinciding with Tin = 0. Considering only 
scalar perturbations from the Einstein-de Sitter Universe, we have at early times 



We can use the residual gauge freedom of the synchronous and comoving gauge to set the 



Starting from these initial conditions, we can assume that the evolution does not 
switch on the off-diagonal components of the PN metric, i.e. that w'^ with a ^ (3 vanish 
at any time. This assumption derives from the physical picture of our one-dimensional 
dynamics. The mass distribution whose self-gravity generates a one-dimensional potential 
is made of parallel sheets of matter. For every point q the peculiar velocity has the same 
direction as the spatial derivative of the peculiar potential, thus the matter moves only 
in the direction perpendicular to the sheets. The collapse of this structure cannot involve 
tensor perturbations, which would lead to the emission of gravitational waves, because it 
cannot undergo any alteration of its shape. However, there is surely a scalar trace part in 
the PN metric, arising from our inflationary initial conditions. In addition, because of the 
asymmetry in the qi spatial direction, the function w\ is assumed to differ from w^^ and 
u^l, whereas the latter functions can only be equal. 

Therefore, the PN expansion is performed according to the following ansatz for the 
metric^ 




(4.10) 



Newtonian perturbation to zero [^]. The initial metric perturbation is therefore given by 
a diagonal PN part. Thus, for the PN perturbation at initial time we have 






722 = 1 + -^h 



(4.13) 



c 



733 = 1 + 



(4.14) 



c 



^Ths indices of the perturbation are lowered (raised) with the background metric 
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with initial conditions 

gm = hin = — j(p. (4.15) 
The PN expansion of the momentum constraint and the Raychaudhuri equation gives 

rj'dih = (dih)' (4.16) 
g" + 2h" + + 2h') + 2,'g' = _ + 2h,^) - (g + 2h) 

We also have the Newtonian energy constraint and evolution equation, connecting the PN 
scalar curvature with Newtonian kinematical quantities: 

-V' + i-'^dfh + 2d,vd,h) = - 1 j (4.18) 

^" + 7^' + i (TT7F ^"^^'^ ^ 2divdih) = . (4.19) 

The momentum constraint is an equation for the spatial derivative of the function h. 
Setting m = dih, it reads 

ri'm = m', (4.20) 
with initial condition = (— 10/3)9i(/3. The solution reads 

m = -ySiVJ (l - • (4.21) 

Then, by spatial integration we obtain 

If) 'i 

h = --^ + -THd,^)' + Co{T), (4.22) 

where the homogeneous mode Cq{t) is a time-dependent constant of integration (w.r.t. 
qi). We can use the Newtonian evolution equation and energy constraint to check the 
consistency of this solution: substitutions of ( [1.22 ) in ( 4.18| ) and ( [4.19| ) leads to the identity. 

Note that the initial condition for the function h sets Co(rin) = 0. In addition, the 
function Cq{t) is an additive term in the perturbation h which would modify the back- 
ground dynamics even in the absence of any initial perturbation (i.e. for ip = 0). Therefore, 
for consistency, we set Co(t) = for all times. 

The Raychaudhuri equation becomes an equation for the function g only, whose solu- 
tion reads 

= . ^ ^r'dMdi^? (,^^. 

^ 3 ^ 5 (-6 + T^dl^) ^ t3 (_6 + ^252^) 21 (-6 + T-^dl^) ' ^ ' 

where C\ and C2 are integration constants. Consistency with our initial conditions, 
Eq. (|4.15| ), requires Ci = 0. 
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In conclusion, the PN metric reads 

.2 \2 , /oi^2/ 



1/10 5 2/0 ^2 

722 = 1 + ^ ( -y(^ + — r2(aiv9)2 

733 = 1 + ^ (- + ^^'(^IV')') • (4.24) 
4.2.1 Determination of the integration constant C2 



In the metric (4.24) the initial condition C2{qi) is still undetermined. To fix it, we take 
advantage of a recent result, obtained in Ref. |]l8|. The authors consider the primordial 
non-Gaussianity set at inflationary epochs on super-Hubble scales. At later times, cosmo- 
logical perturbations re-enter the Hubble radius. They show how the information on the 
primordial non-Gaussianity, set on super-Hubble scales, flows into smaller scale using a 
general relativistic computation. Their calculations are performed in the synchronous and 
comoving gauge, assuming a flat Universe with pressure-less matter plus a cosmological 
constant. They show how the primordial non-Gaussianity affects the PN part of the den- 
sity contrast at second order. Once again the use of inflationary initial conditions in our 
case of globally planar dynamics is justified by our ultimate goal of having a hint on what 
happens in the fully three-dimensional dynamics. 

First of all, for a comparison with the result of Ref. |jl^ our fully non-linear PN 
expression for the density contrast 



6 - T^dl^ c2 \^ 35 (6 _ r^dl^f 3 (6 - r^dl^) 

must be expanded up to second order with respect to the peculiar gravitational potential. 
As usual, we split the density contrast into a first and second order part 5 = d^'^^ + (1/2)5^^^ , 
finding 

Actually, in this expression one can only be sure about the order of the terms that do 
not contain C2, since the latter implicitly depends on the initial peculiar gravitational 
potential, as it will be shown. The first-order term, i.e. T'^dfip/G, obviously coincides with 
the result of linear perturbation theory in the synchronous and comoving gauge and it is 
a Newtonian term, as it is well known. The remaining terms are at least of second order. 
They read 



This expression can be compared with Eq. (45) of Ref. ||18[], by specializing the latter to an 
Einstein-de Sitter background model and to globally planar perturbations^. It reads (in 

^For the Einstein-de Sitter background in (45) of Ref. we set: = 1, f {^Om) ~ 1, Ho = 2/to, 

where the subscript "0" denotes the present time, and D+(r) = t^/tq is the linear growing mode solution. 
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c = 1 units) 

S'-'^ = y - «^^) + Y (2 - aNL) rV(5iV) + ^r^dj^f > (4-28) 

where the deviation of the parameter unl from miity measures the strength of the initial 
(i.e. inflationary) non-Gaussianity (see Ref. [18| for more details). Looking at these expres- 
sions, we first note that the second-order Newtonian term, i.e. (l/18)r^(5^(/?)^, is the same 
in both expression, as it should be. For what concerns the form of C2(gi), it should be 
recalled that it is the initial condition of the PN growing mode cx in the solution (|42^ ). 



Thus, we already know that it must be (at least) a second-order term, since the analogous 
first-order term is Newtonian. The next step is to recognize the PN terms in Eq. ( 4.2^ ). 



Although the explicit powers of c are not shown, one knows from dimensional analysis that 
the second order, i.e. oc <^^, PN terms should be oc and contain two spatial derivatives, 
or they should be oc r'* with four spatial derivatives and so on, in order to have the correct 
powers of c and to be dimensionless, second-order quantities. The PN terms in Eq. ( 4.2^ ) 
are then 

aNL)r\d^y.f (4.29) 



aNL)T\{dfip). (4.30) 



9 

and ^ 

9 

Notice that it is precisely the PN terms which bring all the relevant information about 
(quadratic) primordial non-Gaussianity. Note also that PN terms oc with four spatial 
derivatives and PPN terms are absent (indeed, they would appear at third order in per- 
turbation theory). Now, in the expression ( [4.27 ), it is explicitly shown that (C2/30c'^) t"^ 



is a PN term. Therefore, C2 must contain two spatial derivatives. This fact completely 
determines the form of C2: the most general expression that can be constructed is 

C2 = A{diipf + Bifidfip) . (4.31) 



At this point, notice that the PN terms in Eq. (4.27) 



^^r^9iV (4.32) 



90c2 

and 

:rHd!^? (4.33) 



360c2 

are actually third and fourth-order terms, respectively. Substitution of Eq. ( [4. 311 ) in 
Eq. (14271) leads to 



^ 2/ a ^2| ^ 2/o2 



for our second-order PN terms. By comparison with Eq. (|]2D, we find A and B in terms 
oi qnl we finally obtain 

25 

C2 = y [(1 - 4 (aNL - 1)) (di^f + (4 - 4 (aNL - 1)) ^(df^)] • (4.35) 
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The final expression of our PN metric reads 

5 



. \^ 1 



711= (l--^iV) + o 



6 c 



^^gT^ ( (4 (a^i - 1) - 1) {di^r + (4 (OTVL - 1) - 4) ^dtv ) + 



(6-r2a,V)-|c,(6-r^5,V)^} 
722 = 1 + ^ [-f'P + f^r'idi^f] (4.36) 

1 / 10 5 2/n xS 

733 = 1 + ^ — :rf + T^T (oiv?) 



1/ 10 5 x2 



1/10 5 x2 

These expressions for the metric represent the main result of this paper: they pro- 
vide the post-Newtonian extension of the well-known Zel'dovich solution for plane-parallel 
cosmological dynamics in Newtonian gravity. 

4.2.2 Convergence of the perturbation series 

The actual convergence of the perturbative series requires that the PN metric is much 
smaller than the background Newtonian one. To estimate the order of magnitude of the 
different contributions, one should keep in mind that, on sub-Hubble scales, the peculiar 
gravitational potential is suppressed with respect to the matter density contrast by the 
square of the ratio of the proper scale of the perturbation Xproper to the Hubble radius 
rn = cH^^. Indeed, from the cosmological Poisson equation. 



^^(^)5, (4.37) 



which makes it clear that the gravitational potential divided by the square of the speed 
of light can remain small even on scales characterized by a large density contrast (only 
provided \6\ <C (c/Z^^/Aproper)^), which is indeed at the basis of the well-known validity 
of the Newtonian approach to cosmological structure formation. Let us then consider the 
various terms in the metric of Eq. ( 4.24| ). For 711 



0(1) + 0(5) + 0{6^) + O [0(1) + 0{6) + 0(6^) + 0{6^)] , (4.38) 

where the first three terms belong to the Newtonian part. It is clear that the PN terms are 
sub-leading, being a factor of order O [(p/c^) smaller than the Newtonian ones (provided 
\5\ <C ci? "VAproper)- Similarly, for 722 = 733 we find 



0{l) + 0{^)[0{l) + 0{6)], (4.39) 



thus the same result holds. 



5. A post-Newtonian estimation of kinematical back-reaction 

It has been proposed that the observed increase in the expansion rate of the Universe 
could be due to the back-reaction of the of non-linear sub-horizon cosmic structures on 
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the background Universe expansion [19, 2C, 21 1. The issue of the quantitative relevance of 



back-reaction effects is a largely controversial one (see, e.g. Ref. |22] for a review). It is 
well-known that both standard perturbative treatments even at higher than linear order 
(see e.g. the discussion in Ref. ||2^) and the Newtonian approximation (see in this respect 
the analysis of Ref. ||2^ ) are totally inadequate to correctly evaluate the relevance of back- 
reaction terms in the average Einstein's equations (see also Ref. |24]). For these reasons it 



is tempting to evaluate how PN terms can affect back-reaction. We will then examine this 
issue here using the PN expression of the metric just obtained. The importance of using 
the general relativistic Lagrangian description of cosmological perturbations in connection 
with back-reaction has also been recently stressed in Ref. ||25[| . 

A set of effective Friedmann equations that describe the average dynamics of the real 



inhomogeneous Universe containing irrotational dust have been obtained by Buchert [14| 
by smoothing the Einstein equations by a spatial average for a scalar ^ defined as 

{^)v = ^ I ^Vhd^q, (5.1) 

where h is the determinant of the metric hap and Vv is the volume of the coarse-graining 
comoving domain T) 

Vv= [ Vhd^q. (5.2) 
Jv 

By smoothing the scalar Einstein equations Eqs. ( p. 3D and ( |2.6D , the following effective 
Friedmann equations for the average scale factor ax> = {Vt>/Vvo)^^^ are obtained |14] 

= '-.Gg^Jf (5.3) 



where the source can be viewed as a perfect fluid with effective energy density and pressure 
terms given by 

eff . s QV C^(-^^v 

=^^^^"16^-^6^ 
=-16^ + ^8^- 

obeying the continuity equation 

■eff , o^'D ( eff , Pv\ ^ /c- v^ 

where 

((3)/?)^ = — / ^^'^R^fh S'q (5.8) 
is the average spatial curvature and we have introduced the kinematical back-reaction 

Qv = \m-{&)vf)v-1{^^)v. (5.9) 
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Consistency of Eq. ( ^.T] ) with mass conservation, Eq. ( |2.7D , which can be rewritten as 

Qv + 3—Qv = 0, (5.10) 

requires that the kinematical back-reaction and the mean spatial curvature satisfy the 
integrabihty condition 

{alQvY + c'ai (a|,((%)i,)" = 0. (5.11) 

We remark here that such a condition is a genuinely general relativistic effect, which has 
no analogue in Newtonian gravity, since the curvature ^^^R of comoving hypersurfaces 
vanishes identically in the Newtonian limit. Indeed, in the Newtonian limit the variance 
of the expansion rate and the shear combine to give a total derivative, so Qx> reduces to 



a pure boundary term |23], which becomes negligible if the averaging is performed over 
a volume of the order of the Hubble volume, as should be the case if one is interested 
in the back-reaction of inhomogeneities on the background expansion rate. It was shown 
in Ref. that in order for the back-reaction of cosmological inhomogeneities to drive 
acceleration, the kinematical back-reaction should be Qx> > AttG{q)x>- 

There is another important feature of Eq. ( 5.11| ) which is worth pointing out, namely 
that a factor multiplies the mean spatial curvature term, implying that the Newtonian 
kinematical back-reaction QxJ couples with the PN mean spatial curvature {^^^R)^'^ , the 
PN kinematical back-reaction Q^'^ couples with the post-post-Newtonian mean spatial 
curvature {^^^R)^^^ and so on. 

By performing the conformal rescaling of the metric /ia/3) the velocity-gradient tensor 
splits into a FRW and a peculiar term (recall ( |2.1lD ) 

= i (^^ + , (5.12) 

where a is the scale factor of the underlying FRW background and we use here the conformal 
time coordinate. In terms of the peculiar quantities, the kinematical backreaction becomes 



2 

Explicitly, for our one-dimensional Newtonian background we have CT^ = /3, thus Eq. (|5l^ ) 



becomes 

Qp = -^Wl5, (5.14) 

where a is our Einstein-de Sitter scale-factor, ^ = —Tdl(p/{3{1 — r'^dfip/G)) and (. . .)x> is 
our Newtonian average 

{...)v = £- I ...{l-T^dl^/Q)Sq, (5.15) 

where Vp = a? (l — r'^dfip/G) d?q. We then obtain 

4^2/1 r \ 2 



2aV^ / 1 
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where the last equality follows from the fact that {'d)'D = by construction [23|. 

It follows that any analysis of back-reaction based on the Newtonian approximation 
is irrelevant. In addition, first and second order calculations based on the fully relativistic 
perturbation theory are not sufficient for evaluating back-reaction. This means that fully 
non-linear sub-Hubble perturbations must be taken into account. For these reasons, we 
use our non-linear PN metric to estimate the back-reaction term ( 5.13| ) in the PN approx- 
imation. 

In order to calculate the PN kinematical back-reaction Q^^, we consider the kinemat- 
ical quantities related to the PN peculiar velocity-gradient tensor 

n=^l+it^^^- (5-17) 

In full generality, without any assumption about the Newtonian background, the PN ex- 
pansion of the average integrals in ( ^.13 ) leads to 



~ 3c 

1 / 9 1 1 /-JTO «' 1^/1 1 



CL (X o a 

where w is the trace of the PN perturbation w'^ of the metric and with (. . .)p we indicate 
the Newtonian average (PN corrections in the averaging procedure would yield higher-order 
terms). 

Recalling that in our case a'^ = "& /3, equation ( 5.18 ) reduces to 

= 4(M>©, (5.19) 
where we used the commutation rule for w 

{w)'t, - {w')v = {w^)v - {w)v^)v (5.20) 

and the last equality in Eq. (5.19) follows from the previous Newtonian calculation, having 
neglected the last term, which involves {'d)T>- Finally, using our Newtonian average ( 5.15| ) 
and our solution for W2 , we find 

Qr = -^ I d,{d,^fd\. (5.21) 



81V© Jv 

The PN expansion of (|5.8| ), using our expression for the PN spatial curvature 

(1 — r^a^93/6j 

leads to 

eWv'' = ^ dl^dh. (5.23) 
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It is then clear that, in the case of planar dynamics, both the kinematical back-reaction 
and the mean spatial curvature reduce to purely boundary terms ^ and, as a consequence, 
they cannot lead to acceleration. Also interesting is the fact that {^^^R)^^ tx a~^, which 
provides a self-consistent solution of the integrability condition, Eq. ( ^.IID for vanishing 
Qv- 

5.1 Post-Newtonian corrections to the average scale-factor 

An alternative approach to back-reaction is that of computing directly the PN contribu- 
tion to the average scale-factor. Indeed, using the definition av = (Vd/Vdq)^^^ , we can 
calculate the PN expansion of the average scale- factor, written as 

aD = av + iag^ . (5.24) 



For the Newtonian term, we have for the volume 



Vv = j^a'il- ^djif ] d\ , (5.25) 



thus, neglecting the boundary term, we obtain 

av = — , (5.26) 



which is just the Einstein de Sitter scale factor. 
For the PN term we find 



a^^ = lav {{w)v - {win)v) , (5.27) 

D 



where {■ ■ ■)d is the Newtonian average and our solution gives 



63 

w = 



^\^df^ (d^iff + f r2 (4 (ajvL - 1) + 1) (diipr + (4 (ajvL - 1) + 2) ^df^ 



QOcp 



(5.28) 



6 — r'^djif 
Explicitly, we have 

"^'^ = + [(4(a^,-l) + l)(aiv.)V (5.29) 

+ (4(a,vL-l)-4)vpaiV]|rf'g, 
which can be written as 

ar = S^/ ^^.r^idi^fd'^ (5.30) 
up to negligible boundary terms. The average scale-factor then becomes 

av = avl^l + r^{{l + S)vl)v^ , (5.31) 

vi = —rdiif/S being the (Newtonian) peculiar velocity in the direction xi, which indicates 
a negligible PN correction. 

^Note that these quantities can be seen as boundary terms both in Lagrangian and Eulerian coordinates, 



since 



at the Newtonian level - [dF[q\)/dq )dq = (dF{qi{xi,T)/dx)dx, for any function -F(gi) 
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6. Comparison with the Szekeres solution 



In this paper we have considered the evolution of an irrotational and colhsionless fluid 
in General Relativity in the synchronous and comoving gauge. Following the fluid-flow 
approach, ||l^, it is possible to alternatively describe our system in terms of the fluid 
properties of irrotational dust, i.e. mass density, volume-expansion and shear tensor, and 
the electric and magnetic parts of the Weyl tensor. In addition, in the special case of plane- 
parallel dynamics considered here, the magnetic Weyl tensor vanishes identically, thus 



leading to the so-called silent universe case [26, 27, 28 1, described by the set of equations 



Q = -Oq (6.1) 

e = -6$]2 - -Q^ - AnGg (6.2) 
3 

S = - ^eS - E (6.3) 

E = -3ET. -@E- A-kGqT. , (6.4) 

where is the trace of velocity-gradient tensor, S and E are the eigenvalues of the shear 
tensor and of the electric Weyl tensor 



in the directions q2 and q^. In this framework Croudace et al. |29| obtain what they refere 
to as relativistic Zel'dovich solution, a sub-case of the exact solutions by Szekeres, pO| . 
This solution is made more appealing by the recent growing interest on the Szekeres metric 
in the general framework of studying inhomogeneous cosmologies as possible alternatives 



to FRW (see, e.g. Refs. and fS^] and refs. therein). Ref. |l29[| considers the relativistic 
evolution equations of silent universes. In the special case of local planar symmetry, i.e. 
A2 = A3 = 0, where Aj are the eigenvalues of the matrix dadf^cp, the velocity-gradient tensor 



reads |26| 



0^ = - ( 1 J , (6.6) 

which corresponds to the Newtonian Zel'dovich solution. Croudace et el. then compute 
the associated metric from 

0^ = 2^"" V (6.7) 

via the time-time component of Einstein's equations, i.e. the energy constraint, that closes 
the relativistic fluid-flow equations, completely fixing the spatial dependence of the metric. 
They show that the resulting metric coincides with the Szekeres form 

h22 = /l33 = 1 (6.8) 

hu = (diq^) - ait)c{qi)f (6.9) 

with 

diqa) = dUqi) - \c[qi){{q2f + (93)') " a{t)c{qi) . (6.10) 
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This solution leaves the FRW expansion unperturbed in the directions q2 and 93. In 
fact, following Ref. ||2^, Croudace et al. discard the sub-leading PN trace part of the 
initial conditions ( 4.1C| ), thereby considering perturbations in the direction qi only. On the 



contrary, we kept the one-dimensional initial seed '^{qi) in all directions, thereby allowing 
for perturbations in the component /i22 and /i33 of the metric ah initio. As a consequence, 
our solution for /122 and h^^ changes with time, showing a non-linear dependence on the 
gravitational potential ip. Hence our (approximate) solution cannot be recast into the 
Szekeres form. 

7. Conclusions 

In this paper we have obtained a solution of Einstein's field equations describing the non- 
linear cosmological dynamics of irrotational dust in Lagrangian coordinates for the specific 
case of plane-parallel perturbations, which is exact up to first post-Newtonian order. Our 
solution (Eq. ( 4.36| )), which represents the post-Newtonian extension of the well-known 



Zel'dovich solution for plane-parallel cosmological dynamics in Newtonian gravity, can be 
useful in many respects. Here we have seen just one possible application, in the frame of 
averaging and back-reaction of cosmological inhomogeneities. In a forthcoming paper we 
will analyze the properties of photon geodesies in our metric, to study, for instance, the 
luminosity distance-redshift relation in a strongly inhomogeneous Universe such as the one 
described by our metric. 

Our results also allow us to study the final stages of plane-parallel collapse, which ob- 
viously leads to a shell-crossing, pancake-like singularity. Caustic formation is considered 
the main limitation of the Lagrangian description, both in the Newtonian approximation 
and in General Relativity. As well-known, caustics arise because several fluid elements 
coming from different positions may converge to the same Eulerian position, thus forming 
infinite density regions. Such a pathological behavior occurs in our case when / = — 1, 
i.e. when r'^dlip/Q = 1: at this time the determinants of both the Newtonian ( [4. 61) and 
PN metric (g^) go to zero, while the density contrast at both the Newtonian, Eq. (fl.^), 
and PN level, (|^), becomes infinite, as does the PN spatial curvature in Eq. ( ^.22] ). The 



appearance of shell-crossing singularities can be understood as indicating the breakdown 
of the dust approximation, rather than the occurrence of a true physical singularity of the 
gravitational collapse. The formation of caustics appears as an artifact of the extrapolation 
of the pressure-less fluid approximation beyond the point at which pressure has become 
important. An important result of our analysis is that the divergence of the PN spatial cur- 
vature at caustic formation is completely eliminated by the spatial smoothing procedure. 
Indeed TZ^^ diverges like {1 + 5), which is exactly compensated by the square root of the 
spatial metric determinant oc 1/(1 + 6). This very fact confirms that the instability found 
in Ref. using a gradient expansion technique, and in Ref. [^, using a different ap- 
proximation scheme, cannot be interpreted as a consequence of a shell-crossing singularity, 
but really arises from the back-reaction of sub-Hubble modes. 

Our results here bring both good and bad news for the back-reaction of cosmic inho- 
mogeneities to represent a potentially viable alternative to dark energy at a fundamental 
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level. Let us start with the bad news: as we have shown the kinematical back-reaction 
scalar remains a negligible boundary term also at the post-Newtonian level, so that no rel- 
evant back-reaction effect is implied by our PN solution. We should stress, however, that 
such a result is clearly a consequence of our very specific (but analytically solvable) model, 
which relies on one-dimensional dynamics. It is a very reasonable hypothesis that such a 
result will be modified by a (more complex!) full 3D calculation; this will be the subject 
of a subsequent paper. The good news are: the very fact that the Lagrangian approach 
allows to obtain a quantitative estimate of back-reaction and the fact that the appear- 
ance of shell-crossing singularities does not lead to divergences in the average dynamics of 
inhomogeneous dust Universes. 

Also interesting, in our opinion, is the fact that the PN metric depends on the pri- 
mordial non-Gaussianity strength parameter cnl (which is related to the best known non- 
Gaussianity parameter /nl by /nl = (5/3)aNL! for large onl)- Although in our final results 
for 1 {^^'^R)^^ and this quantity disappeared (in that it only enters in boundary 
terms), one can safely predict that this parameter will affect back-reaction quantities for 
general three-dimensional perturbations, hence suggesting the intriguing possibility of a 
potential quantitative link between primordial non-Gaussianity and back-reaction. 
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Abstract: We study the general relativistic non- linear dynamics of self-gravitating ir- 
rotational dust in a cosmological setting, adopting the comoving and synchronous gauge, 
where all the equations can be written in terms of the metric tensor of spatial hyper-surfaces 
orthogonal to the fluid flow. Performing an expansion in inverse powers of the speed of 
light, we obtain the post-Newtonian equations, which yield the lowest-order relativistic 
effects arising during the non-linear evolution. We then specialize our analysis to globally 
plane-parallel configurations, i.e. to the case where the initial perturbation field depends 
on a single coordinate. The leading order of our expansion, corresponding to the "Newto- 
nian background" , is the Zel'dovich approximation, which, for plane-parallel perturbations 
in the Newtonian limit, represents an exact solution. This allows us to find the exact 
analytical form for the post-Newtonian metric, thereby providing the post-Newtonian ex- 
tension of the Zel'dovich solution: this accounts for some relativistic effects, such as the 
non-Gaussianity of primordial perturbations. An application of our solution in the context 
of the back-reaction proposal is eventually given, providing a post-Newtonian estimation 
of kinematical back-reaction, mean spatial curvature, average scale-factor and expansion 
rate. 
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1. Introduction 

Large-scale structures in the Universe grew by gravitational instability around primordial 
seed perturbations generated during inflation. In this context, it is assumed that the 
matter mainly consists of Cold Dark Matter (CDM), which, prior to caustic formation, 
behaves like a perfect fluid of dust. The study of the gravitational dynamics of CDM can 
be performed by different techniques and approximations, depending on the specific range 
of scales under investigation. 

The most widely used approximation is standard perturbation theory (see e.g. Ref. Q), 
which is an expansion in powers of the amphtude of fluctuations around a Friedmann- 
Robertson- Walker (FRW) background. The range of applicability of this perturbation 
technique is that of small fluctuations around a FRW background, with no extra limitations 
on scale. Going to higher perturbative orders generally helps to follow the gravitational 
instability process on a longer time-scale and to account for non-linear phenomena, i.e. to 
resolve smaller scales. 

Another important approximation scheme is the Newtonian one, which is known to 
produce accurate results as long the gravitational fleld is weak and the motion of both 
sources (i.e. CDM fluid elements) and test particles is slow. In the cosmological framework. 
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the range of validity of the Newtonian approximation concerns scales much larger than the 
Schwarzschild radius of collapsed bodies and much smaller than the Hubble horizon (see 
e.g. Ref. 0). This approximation is usually carried out in the Poisson gauge, where 
the FRW metric is perturbed by a single lapse function which represents the peculiar 
gravitational potential sourced by matter density fluctuations. This approach provides 
the well-known Newtonian Eulerian description of cosmological structure formation in the 
expanding universe. The extension of this approximation to deal with fast moving test 
particles (like photons) is the so-called weak-field limit, where a scalar potential is accounted 
for in the spatial part of the line-element. 

In the context of Newtonian gravity, departures from the FRW background may be also 
analyzed in the Lagrangian frame. In the unperturbed situation, observers comoving with 
the Hubble flow are Lagrangian, but when density fluctuations are taken into account, they 
become Eulerian observers with respect to the peculiar motion. The equations describing 
the dynamics of CDM, which we model here through irrotational dust, are written in La- 
grangian coordinates, in terms of the displacement vector S from Lagrangian to Eulerian 
coordinates. As in the Eulerian case, it is impossible to work out the general analytical 
solution for S and a perturbative approach is again introduced by means of an expansion 
in powers of the displacement vector, the background being once more represented by the 
FRW model. The linear result is the so-called Zel'dovich approximation, Q. The peculiar- 
ity of this treatment, at any order, is that, while the displacement vector is calculated from 
the equations at the required perturbative order, all the other dynamical variables, such as 
the mass density, are calculated exactly from their non-perturbative definition. Since the 
equations in Lagrangian coordinates are intrinsically non-linear in the density, what comes 
out is a fully non-linear description of the system, which, though not being generally exact, 
"mimics" the true non-linear behavior. This perturbation treatment basically exploits the 
advantages of the Lagrangian picture, leading, in particular, to a more accurate descrip- 
tion of high density regions. Its limitations are generally set by the emerging of caustic 
singularities. 

The Newtonian analysis can be improved by a post-Newtonian (FN) approach. In the 
Newtonian approximation two conditions must be satisfied: weak gravitational field and 
slow motion. Relaxing one of these two conditions leads to two different improvements 
of the Newtonian description. The FN approximation is suitable for a system of slowly 
moving particles bound together by gravitational forces and can be useful in cosmology 
(see Refs. gi ), as it can be used to account for the moderately non-linear gravitational 
field generated during the highly non-linear stage of the evolution of matter fiuctuations. 
The FN approximation in the cosmological framework has been studied in Eulerian (or 
more generally non-comoving) coordinates in Refs. [^, ^, ^ |^, while in Ref. [11 1 a hy- 
brid approximation scheme is proposed which upgrades the weak-field limit of Einstein's 
field equations to account for post-Newtonian scalar and vector metric perturbations and 
for leading-order source terms of gravitational waves, while including also the first and 
second-order perturbative approximations. An alternative approach has been recently pro- 
posed in Ref. [12|, whose authors introduce an effective fiuid description for small-scale 
non-linearities in the framework of General Relativity. The FN approximation to cosmo- 
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logical structure formation within the Lagrangian picture has been studied by Matarrese 



& Terranova in Ref. |13], whose approach will be the basis for our analysis here. 

The aim of this paper is to find a PN metric describing the non-linear stages of the 
gravitational instability in an Einstein-de Sitter universe, adopting the Lagrangian picture 
i.e. the synchronous and comoving gauge. In this gauge the PN approximation is formally 
obtained via an expansion in inverse powers of the speed of light c. We choose the lead- 
ing order solution of our expansion to be the Newtonian metric related to the Zel'dovich 
approximation, in the very special case of globally planar dynamics due to a perturbation 
seed (/9 depending on a single Lagrangian coordinate. 



Our approach here is very different from those of Refs. [ITl, 15|, who proposed rela- 



tivistic generalizations of the Zel'dovich approximation. Ref. |14| introduced a relativistic 



tetrad-based, perturbative approach, which is then solved to linear order and used to ob- 
tain non-perturbative expressions for the velocity-gradient tensor and mass density. The 



solution of Ref. ||15|| is instead equivalent to a relativistic second-order perturbation theory 
treatment in the synchronous and comoving gauge, in which all quantities (metric, velocity- 
gradient tensor and mass density) are calculated at second order, thereby partially missing 
the non-perturbative character of the Zel'dovich approximation. Our approach instead aims 
at obtaining a non-perturbative description of both metric and fluid properties (velocity- 
gradient tensor and mass density), within the post-Newtonian approximation of General 
Relativity: our expansion in inverse powers of the speed of light is fully non-perturbative 
from the point of view of standard perturbation theory; thus our results contain all sec- 
ond and higher-order terms of standard perturbation theory calculations, as long as they 
are post-Newtonian and one deals with the plane-parallel dynamics. In our approxima- 
tion scheme the Zel'dovich solution represents the Newtonian background over which PN 
corrections can be computed as small perturbations. 

The plan of the paper is as follows. In Section 2 we introduce the relativistic equations 
governing the dynamics of irrotational dust. Section 3 deals with the Newtonian limit. 
Section 4 is devoted to the PN approximation: we first describe the Zel'dovich Newtonian 
background chosen as a basis for our PN expansion, showing explicitly all the dynamical 
quantities involved. Then, we obtain and solve analytically the PN equations for the metric. 
In particular, we show that our solution accounts for the non-Gaussianity of primordial 
cosmological perturbations. To conclude this section, we briefly discuss the accuracy of 
our perturbative expansion. Section 5 deals with an application of our solution: we use 
our metric in the context of the back-reaction proposal. To be more specific, we give an 
estimate of the PN contribution to the kinematical back-reaction, mean spatial curvature. 



average scale-factor and expansion rate, as defined in Refs. |16, 17]. Section 6 provides a 
comparison between our PN solution and the Szekeres metric. Conclusions are drawn in 
Section 7. 

2. Relativistic dynamics of irrotational dust in the Lagrangian picture 

In this section we recall the relativistic equations governing the evolution of irrotational 
dust following the Arnowitt-Deser-Misner (ADM) formalism, [|^]. Calculations are made 
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in the synchronous and comoving gauge. In this gauge both the world-lines of the observers 
coincide with the geodesies of the matter and the proper time along the geodesies coincides 
with the cosmic time. This reference frame is therefore Lagrangian. Actually, the possibility 
of making these two gauge choices simultaneously is a peculiarity of irrotational dust, which 
holds also in the non-linear stage. 
The line element is given by ^ 

ds^ = -(?dt^ + Kp{t, (\)dq'^dq^ , (2.1) 

where Qa are the spatial Lagrangian coordinates and t is the proper time of our fluid 
elements. The matter stress-energy tensor reads T"^ = q(?u'^v!' with q the mass density 
and the fluid four-velocity (normalized to u^Ua = —1), with components = (1, 0, 0, 0). 
The ADM equations can be written in terms of the velocity-gradient tensor defined as 

= cu^ = \k^-Kp , (2.2) 

where dots denote partial differentiation with respect to the time t and semi-colons stand 
for covariant differentiation. The tensor 0^ represents the extrinsic curvature of constant 
t hyper-surfaces. 

The 10 Einstein equations split into 4 constraints and 6 evolution equations. The 
time-time component of the Einstein equations is the so-called energy constraint of the 
ADM approach, which reads 

- e^G^ + = IQttGq , (2.3) 

where the volume-expansion scalar is the trace of the velocity-gradient tensor and ^^^R is 
the trace of the 3-dimensional Ricci tensor ^^^R^ of the spatial hyper-surfaces of constant 
time. 

The space-time components give the momentum constraint, 

Q%a = ■ (2-4) 

The space-space components govern the evolution of the velocity-gradient tensor and rep- 
resent the only truly evolution equations. They read 

8^ + 69^ + = 47rGe5^ . (2.5) 

Taking the trace of the latter equation and combining it with the energy constraint yields 
the Raychaudhuri equation for the volume- expansion scalar 0, 

e + 2S2 + le2 + 47rGe = 0, (2.6) 

where = G°— (1/3)B(^" is the s/iear tensor and T? = (1/2)S^Sq its magnitude. Finally, 
the continuity equation for the stress-energy tensor, T^p = QUaUp, leads to 

Q = -Qq. (2.7) 

^We use Greek letters for purely spatial indices, a, p... = 1, 2, 3, whereas Latin ones a, h.. = 0, 1, 2, 3 
label space-time indices. 
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whose solution reads 



h being the determinant of the spatial metric. Here and in what follows quantities with a 
subscript "in" are meant to be evaluated at some initial time tjn. 

2.1 Conformal rescaling and FRW background subtraction 

With the purpose of studying gravitational instability in the Einstcin-de Sitter background, 
we first factor out the homogeneous and isotropic expansion of the Universe. To this aim, 
we perform a conformal rescaling of the metric with the conformal factor a(r) oc r^, the 
scale factor of the Einstein de-Sitter model. This procedure involves the transformation of 
the proper time t to conformal time r via dr = dt/a(t). The line-element is written in the 
form 

ds'^ = a^{T) \-c^dT^ + 7a/3(r,q)(ig°dg^] , (2.9) 

where 7^/3 (t, q) = h^g (t (r) , q) /a^ (t (r)) is the conformal spatial metric. In order to factor 
out the Einstcin-de Sitter background, the isotropic Hubble flow, with velocity-gradient 
tensor (a'/a)(5^, is subtracted from the velocity-gradient tensor, leading to a tensor which 
describes gradients of the peculiar velocity only: 

^« = acu'^ - ^<5^ = ^j^H^ , (2.10) 

where n° = (1/a, 0, 0, 0) and primes denote partial differentiation with respect to r. The 
tensor represents the extrinsic curvature of constant r spatial hypersurfaces. The 

matter content is also written in terms of the density contrast, defined as the dimen- 
sionless deviation of the matter density from that of the Einstein-de Sitter background, 
ft(T) = 3/(2irGaV2), 

i(r.q) (2.11) 

Qb{T) 

The ADM equations are then recast in a more convenient form describing the evolution of 

the peculiar velocity-gradient tensor. The energy and momentum constraints become 

o 24 
'd'^ -■d'^'d^ + -'d + c^n = —S (2.12) 

= , (2.13) 

where Va is the covariant derivative associated with 

is the conformal three-dimensional (scalar) curvature. The evolution equation becomes 

7?^' + ^7?^ + + ^7?<5^ + c^n'l = ^66^ . (2.14) 
After replacing the density from the energy constraint, this equation is written as 

+^[47^^-7^<5^]. (2.15) 
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The Raychaudhuri equation describes the evolution of the pecuhar volume expansion scalar 
and reads 

^' + _^9 + ^/;^^ + _j = 0. (2.16) 



Finally, the solution of the continuity equation ( |2.8| ) is written in terms of the contrast 6 
as 

5(q,t) = (i + Mq))^^^-i, (2.17) 

where 7 is the determinant of the conformal spatial metric 7^/3. 

The main advantage of this formalism is that there is only one dimensionless variable in 
the equations, namely the spatial metric tensor 7^^, which is present with its partial time 
derivatives through and with its spatial gradients through the spatial Ricci tensor 7^^. 
A relevant advantage of having a single tensorial variable, for the following PN expansion, 
is that there can be no extra powers of c hidden in the definition of different quantities. 

3. Newtonian approximation 

The Newtonian equations in Lagrangian form can be obtained from the full relativistic 
equations of the previous section by an expansion in inverse powers of the speed of light. 
As a consequence of our gauge choice, however, no odd powers of c appear in the equations, 
so the expansion parameter is 1/c^. The spatial metric is then expanded in the form 

laP=lal3 + -^Wal3 + 0(^^^ . (3.1) 

Let us first concentrate on the Newtonian metric Ja/s- The Newtonian limit (c — )• oo) of 
the energy constraint and evolution equation requires that the spatial Ricci tensor is zero. 
Thus, the Newtonian result of vanishing spatial curvature is recovered (see Refs. |19, 0]). 
This important conclusion implies that 7q,^ can be transformed to the Euclidean metric 
6ai3 globally. In other words, at each time r there exist global Eulerian observers comoving 
with the Hubble flow for which the components of the metric are Sa^. This means that, 
according to the tensor transformation law, we can write 

7q/3 = Sf^uJ^Jp , (3.2) 

where 

is the Jacobian of the transformation 

x(q,r) =q + 5(q,T) (3.4) 

and x^, = 1,2,3), are the Eulerian (comoving with the isotropic FRW background) co- 
ordinates. The tensor dS^/dq°' is called "deformation tensor". All the information about 
the motion is contained in the mapping between the Eulerian and Lagrangian coordinates, 
whose evolution is described by the Raychaudhuri equation and the momentum constraint. 
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In fact, contrary to the evolution equation and the energy constraint, the Raychaudhuri 
equation and the momentum constraint contain no exphcit powers of c, thus preserve their 
form in going to the Newtonian Hmit. These equations determine the background Newto- 
nian metric 7q,/3: they are indeed the Newtonian equations in the Lagrangian approach to 
irrotational dust dynamics (see Ref. [|l9|). Thus, it makes sense to recast these equations in 
terms of the map between comoving Eulerian and Lagrangian coordinates. One can equiv- 



alently write the equations in terms of either the displacement vector, as in Refs. [2C, 21 



and deformation tensor or in terms of the Jacobian matrix of the map, as in Ref. |13|. 



Here, following Ref. [13], we consider the Newtonian equations in terms of the Jacobian 



matrix ([3.31). They read 



= (3.5) 



e''«V;j;.-, = 0, (3.6) 

where J is the determinant of the Jacobian matrix. We assumed for simplicity = and 
used the residual gauge freedom of the synchronous and comoving gauge to set = 1 in 
the Newtonian limit, as in Ref. | p!3[ . 

It is impossible to follow the dynamics in an exact analytical way and one is forced to 
use approximation techniques to solve these equations. Here the perturbative expansion is 
performed in terms of the Jacobian matrix. The key point is that a slight perturbation of 
the Lagrangian particle paths carries information on the non-linear dynamics, as may be 
seen directly from the expression of the density contrast 

5 = 1-1, (3.7) 

which shows that, even if the displacement vector is small, the corresponding density con- 
trast could be large. In other words, already the first-order Lagrangian equations describe 
the mildly non-linear regime of the gravitational instability. 

Expanding and solving the previous equations to first order, the well-known Zel'dovich 
approximation is recovered: 

x(q,r) = q + D(r)V<I>(q), (3.8) 

where -D(t) oc a(T) oc is the growing mode solution for the Einstein-de Sitter model and 
$(q) has to be ascribed to the chosen initial conditions. 



4. Post-Newtonian approximation 

4.1 Characterization of the Newtonian background 

The starting point of our PN expansion is the Newtonian background described by the 
Zel'dovich approximation, with the peculiar gravitational potential depending on the con- 
formal time and on the Lagrangian coordinate q\ only. As it is well known, in the particular 
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case of planar perturbations the Zel'dovich approximation yields an exact solution of the 
Newtonian equations. The Zel'dovich solution is given by 

xi = qi + D{T)di'^in{qi) 
X2 = q2 

X3 = q3- (4.1) 

where di = d/dqi and the potential <I>in is defined so that 5f $in = — ^in/An and is related to 
the initial peculiar gravitational potential (p by the cosmological Poisson equation, yielding 

^■^- = -T7^2 • (4-2) 



The Jacobian matrix is 



1 - T^^iV/e 

Jp = \ 1 I • (4.3) 

For the following calculations, it is useful to define the function / 

/ := D{T)dl^,n = -^dfif (4.4) 
b 

and the function r] = ln(l + /). Hereafter, the peculiar gravitational potential is meant 
to be evaluated at the initial time Tin and the subscript "in" is dropped for notational 
convenience. 

The components of the metric in Lagrangian coordinates are given by (recall Eq. ( |3.2| )) 

7a/3 = Sau. K + D{T)d''da'^] [5^ + D{T)d^dp^] , (4.5) 

where we used the fact that at first order in the displacement vector covariant and partial 
derivatives with respect to the coordinates qa coincide, since the Christoffel symbols are 
second-order quantities. In our case $in = $in(gi) and for the Zel'dovich metric we find 

7a/3 = I 1 1 (4.6) 

or in more compact form 

/(! + /)' \ 
7a/3 = I 1 1 . (4.7) 

It is important to keep in mind that this Newtonian metric is non-linear with respect 
to the peculiar gravitational potential, thus it characterizes the mildly non-linear stage of 
the gravitational instability. Starting from this metric at first order in the displacement 
vector, all the other dynamical variables are calculated exactly. The only non-vanishing 
component of the peculiar velocity-gradient tensor is "di = tj' and for the shear tensor we 
have a\ = 2ri'/3 and = al = —a\/2. 
Finally the density contrast ( |3.7] ) takes the form 

1 
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4.2 Post-Newtonian expansion 

For the PN expansion, we write the metric in the form 

laf} = lap + • (4-9) 

For consistency with our Newtonian background solution, in which the pecuhar gravita- 
tional potential depends only on qi, the PN perturbation Wap can be assumed to depend 
on the conformal time and on the Lagrangian coordinate qi only. 

A clue for the form of the perturbation w^p follows from the initial conditions of 
cosmological perturbations. Even though our case of globally planar dynamics is purely a 
toy-model we prefer to assume that our initial perturbation seed is consistent with having 
been generated during inflation in the early Universe, so that our analytical results can 
give us a hint of what happens in the real Universe, where the initial perturbation seed 
is a random field that depends on all the spatial coordinates. We then set our initial 
conditions at the end of inflation, effectively coinciding with Tin = 0. Considering only 
scalar perturbations from the Einstein-de Sitter Universe, we have at early times 

7a/3 = ~ ^'^) "^"^ ~ \dc,dpip . (4.10) 

We can use the residual gauge freedom of the synchronous and comoving gauge to set the 
Newtonian perturbation to zero (see Ref. |l^). The initial metric perturbation is therefore 
given by a diagonal PN part. Thus, for the PN perturbation at initial time we have 

10 
3 

Starting from these initial conditions, we can assume that the evolution does not 
switch on the off-diagonal components of the PN metric, i.e. that w'^ with a ^ (3 vanish 
at any time. This assumption derives from the physical picture of our one-dimensional 
dynamics. The mass distribution whose self-gravity generates a one-dimensional potential 
is made of parallel sheets of matter. For every point q the peculiar velocity has the same 
direction as the spatial derivative of the peculiar potential, thus the matter moves only 
in the direction perpendicular to the sheets. The collapse of this structure cannot involve 
tensor perturbations, which would lead to the emission of gravitational waves, because it 
cannot undergo any alteration of its shape. However, there is surely a scalar trace part in 
the PN metric, arising from our inflationary initial conditions. In addition, because of the 
asymmetry in the qi spatial direction, the function w\ is assumed to differ from w^^ and 
u^l, whereas the latter functions can only be equal. 

Therefore, the PN expansion is performed according to the following ansatz for the 
metric^ 

711 = {l + f? + \{l + f?g (4.12) 
722 = 1 + \h (4.13) 
733 = 1 + ^/1, (4.14) 



^/3in = --^^0 ■ (4.11) 



^The indices of the perturbation are lowered (raised) with the background metric (7° 
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with initial conditions 

9in = hn = — -^f- (4.15) 

The PN expansion of the momentum constraint and of the Raychaudhuri equation 
gives 

r^'dih = {dih)' (4.16) 

r 1 + / 

We also have the energy constraint and evolution equation, connecting the PN scalar 
curvature with Newtonian kinematical quantities: 

+ {-2dlh + 2d,^d,h) = I - l) (4.18) 

^" + 7^' + \ {llfY ("^^1^ + 2diiidih) = . (4.19) 

The momentum constraint is an equation for the spatial derivative of the function h. 
Setting m = dih, \t reads 

r^'m = m! , (4.20) 
with initial condition = (— 10/3)9i99. The solution reads 

m = -^di^ {l - ^SiV) • (4.21) 

Then, by spatial integration we obtain 

in c; 

h = --^ + -T\d,^f + C,{T), (4.22) 

where the homogeneous mode Co(t) is a time-dependent constant of integration (w.r.t. 
qi). We can use the Newtonian evolution equation and energy constraint to check the 
consistency of this solution: substitutions of ( [1.22D in ( |4.18| ) and ( [4.19| ) leads to the identity. 

Note that the initial condition for the function h sets Co (Tin) = 0. In addition, the 
function Co(t) is an additive term in the perturbation h which would modify the back- 
ground dynamics even in the absence of any initial perturbation (i.e. for = 0). Therefore, 
for consistency, we set Co(r) = for all times. 

The Raychaudhuri equation becomes an equation for the function g only, whose solu- 
tion reads 

= T^g2 , Ci br^dl^jd,^? 

^ 3 ^ 5 (-6 + T^dl^) ^ (-6 + T^djip) 21 (-6 + T^d^ip) ' ^ ' 

where Ci and C2 are integration constants. Consistency with our initial conditions, 
Eq. (|4.15| ), requires Ci = 0. 
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In conclusion, the PN metric reads 

/ t2 \2 f2lT^C2-25T''dMdw?-S50ip{-6 + T''dl^)' 

711 = ( 1 - — ) + — (-6 + T ' 



6 ' W ' '"M 3780 



1/10 5 2/n x2 

722 = l + ^i- — V+—T\diipf 



733 = 1 + ^ ( -y ^ + ^r\d,^f ] . (4.24) 



4.2.1 Determination of the integration constant C2 



In the metric ( 4.24 ) the initial condition 6*2(^1) is still undetermined. To determine it we 



take advantage of the results obtained in Refs. [22, 23]. The authors consider the primor- 
dial non-Gaussianity set at inflationary epochs on super-Hubble scales. At later times, 
cosmological perturbations re-enter the Hubble radius. They show how the information 
on the primordial non-Gaussianity, set on super-Hubble scales, flows into smaller scale 
using a general relativistic computation. Their calculations, which are performed in the 
synchronous and comoving gauge, show how the primordial non-Gaussianity affects the 
PN part of the density contrast at second order. Once again the use of inflationary initial 
conditions in our case of globally planar dynamics is justified by our ultimate goal of having 
a hint on what happens in the fully three-dimensional dynamics. 

First of all, we consider our fully non- linear PN expression for the density contrast. 
The PN contribution is given by (5^^ = (1/2) (l + 6) {w-m — w), where w is the trace of the 
PN perturbation w'^ of the metric 

_ -315O99 {r'^dlf - 6) + t2 [63C2 + 50{di<ff {ir'^dfip - 2l)] 

"'^ 315 (r^a^V-e) • 



For a comparison with the result of Ref. |^3[| , our expression for the density contrast 

6 - T^dl^ ^ 35 (6 - T^dl^f 3 (6 - r^dlip) ) ^ " 

must be expanded up to second order with respect to the peculiar gravitational potential. 
As usual, we split the density contrast into a first and second order part 5 = d^'^^ + (1/2)5^^^ , 
finding 

Actually, in this expression one can only be sure about the order of the terms that do 
not contain C2, since the latter implicitly depends on the initial peculiar gravitational 
potential, as it will be shown. The first-order term, i.e. r^S^i/j/G, obviously coincides with 
the result of linear perturbation theory in the synchronous and comoving gauge and it is 
a Newtonian term, as it is well known. The remaining terms are at least of second order. 
They read 



30?" +90?" ^1^-9?" (^^^^ + 360?" (^^^^ +I8" (^^^^ 
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This expression can be compared with Eq. (45) of Ref. ||23[], by speciahzing the latter to an 
Einstein-de Sitter background model and to globally planar perturbations^. It reads (in 
c = 1 units) 

S^'^ = y (I - -nl) r\d,^? + y (2 - a^,) r^dh) + Y^Adl^f ^ (4-29) 
where the deviation of the parameter unl from unity measures the strength of the initial 



(i.e. inflationary) non-Gaussianity (see Ref. [23| for more details). Looking at these expres- 
sions, we first note that the second-order Newtonian term, i.e. (l/18)r^(5^(/?)^, is the same 
in both expression, as it should be. For what concerns the form of C2{qi), it should be 
recalled that it is the initial condition of the PN growing mode oc in the solution (|42^ ). 



Thus, we already know that it must be (at least) a second-order term, since the analogous 
first-order term is Newtonian. The next step is to recognize the PN terms in Eq. (l4.2g| ). 



Although the explicit powers of c are not shown, one knows from dimensional analysis that 
the second order, i.e. oc ip^, PN terms should be oc and contain two spatial derivatives, 
or they should be oc with four spatial derivatives and so on, in order to have the correct 
powers of c and to be dimensionless, second-order quantities. The PN terms in Eq. ( 4.29] ) 
are then 

^i{^^-^NL)r\d,^f (4.30) 

and 

^(2-a^i)rV(5iV). (4.31) 

Notice that it is precisely the PN terms which bring all the relevant information about 
(quadratic) primordial non-Gaussianity. Note also that PN terms oc with four spatial 
derivatives and PPN terms are absent (indeed, they would appear at third order in per- 
turbation theory). Now, in the expression ( [4.28| ), it is explicitly shown that (C2/30c'^) t"^ 
is a PN term. Therefore, C2 must contain two spatial derivatives. This fact completely 
determines the form of C2: the most general expression that can be constructed is 

C2 = A{di^f + B^{dl^) . (4.32) 



At this point, notice that the PN terms in Eq. (4.28) 



'^2 ,4 02 



T^di^ (4.33) 



90c2 

and 

' ^r^dl^f (4.34) 



360c2 

are actually third and fourth-order terms, respectively. Substitution of Eq. ( [4. 321 ) in 



Eq. (4.28) leads to 



^ ^ ^ 2/0 \2 , ^ 2 /a2 



T\d,^Y + — -rV(5iV) (4.35) 



30c2 9c2/ ' ^' 30c2 



^For the Einstein-de Sitter background in (45) of Ref. we set: flom = 1, /(f^Om) ~ 1, 'Ho = 2/to, 
where the subscript "0" denotes the present time, and D+{t) = /tq is the linear growing mode solution. 
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for our second-order PN terms. By comparison with Eq. ( [4.29 ), we find A and B in terms 
of otvl we finally obtain 
25 

C2 = y [(1 - 4 (ajvL - 1)) [diiff + (4 - 4 {aNL - 1)) ^{djif)] . (4.36) 
The final expression of our PN metric reads 

711 = (1 - y 9iV) + ^ { -^r^ ((4 (a^vL - 1) - 1) (di^? + (4 (a^vL - 1) - 4) ipdfip) + 



(6-r2a,V)-|^(6-r^a,V)^} 



722 = 1 + ^ (-y + ^r2(ai<^)' ) (4.37) 



1/10 5 2/n x2 

733 = 1 + ^ ( -y¥'+ (9lV9) 

These expressions for the metric represent the main result of this paper: they pro- 
vide the post-Newtonian extension of the well-known Zel'dovich solution for plane-parallel 
cosmological dynamics in Newtonian gravity. 

4.2.2 Convergence of the perturbation series 

The actual convergence of the perturbative series requires that the PN metric is much 
smaller than the background Newtonian one. To estimate the order of magnitude of the 
different contributions, one should keep in mind that, on sub-Hubble scales, the peculiar 
gravitational potential is suppressed with respect to the matter density contrast by the 
square of the ratio of the proper scale of the perturbation Xproper to the Hubble radius 
rff = cH~^. Indeed, from the cosmological Poisson equation, 

^-(^p, (4.38) 



which makes it clear that the gravitational potential divided by the square of the speed 
of light can remain small even on scales characterized by a large density contrast (only 
provided \6\ <^ (c-ff^^/Apropcr)^)i which is indeed at the basis of the well-known validity of 
the Newtonian approach to cosmological structure formation. We must also recall that, in 
the Newtonian limit, the square of the peculiar velocity is of the same order as the peculiar 
gravitational potential and both remain small even in the non-linear regime of structure 
formation. 

Let us then consider the various terms in the metric of Eq. ( [4.24| ). For 711 we find 

O ( +o(^)+o( , (4.39) 

where the first term belongs to the Newtonian part. Similarly, for 722 = 733 we have 

0(l) + 0(^). (4.40) 
It is clear that the PN terms are sub-leading. 
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5. A post-Newtonian estimation of kinematical back-reaction 



It has been proposed that the observed increase in the expansion rate of the Universe 
could be due to the back-reaction of the of non-hnear sub-horizon cosmic structures on 



the background Universe expansion p4|, 25, B6|. The issue of the quantitative relevance of 



back-reaction effects is a largely controversial one (see, e.g. Ref. |27] for a review). It is 



well-known that both standard perturbative treatments even at higher than linear order 



(see e.g. the discussion in Ref. |25|) and the Newtonian approximation (see in this respect 
the analysis of Ref. [|2^ ) are totally inadequate to correctly evaluate the relevance of back- 
reaction terms in the average Einstein's equations (see also Ref. ||2^). For these reasons it 
is tempting to evaluate how PN terms can affect back-reaction. We will then examine this 
issue here using the PN expression of the metric just obtained. The importance of using 
the general relativistic Lagrangian description of cosmological perturbations in connection 
with back-reaction has also been recently stressed in Ref. pO[ |. 

A set of effective Friedmann equations that describe the average dynamics of the real 
inhomogeneous Universe containing irrotational dust have been obtained by Buchert in 
Ref. [^] by smoothing the Einstein equations by a spatial average for a scalar ^ defined 
as 

mv = ^ / ^Vhd\, (5.1) 

where h is the determinant of the metric ha/s and V-p is the volume of the coarse-graining 
comoving domain D 

Vv= [ Vhd^q. (5.2) 
Jv 



By smoothing the scalar Einstein equations, Eqs. (2^) and (|2.6D , the following effective 
Friedmann equations for the average scale factor a© = (V-p/Vc,,)^''^ are obtained 



-kGq^' (5.3) 

where the source can be viewed as a perfect fluid with effective energy density and pressure 
terms given by 



V 

obeying the continuity equation 



=-16^ + ^8^- 



where 



^?/' + 3?U^^^ + 4^)=0, (5.7) 



((%)^ = 7^ / ^^^RVh d^q (5.8) 
Vv Jv 
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is the average spatial curvature and we have introduced the kinematical back-reaction 

Qv = l{{e-{e)vf)v-2{^^)v. (5.9) 



Consistency of Eq. ( ^.7| ) with mass conservation, Eq. (|2.7|), which can be re-written as 

Qj, + 3—Qv = 0, (5.10) 
a© 

requires that the kinematical back-reaction and the mean spatial curvature satisfy the 
integrability condition 

{a%Qvy + c'ai (a|,((%)i,)" = 0. (5.11) 

We remark here that such a condition is a genuinely general relativistic effect, which has 
no analogue in Newtonian gravity, since the curvature ^^^R of comoving hypersurfaces 
vanishes identically in the Newtonian limit. Indeed, in the Newtonian limit the variance 
of the expansion rate and the shear combine to give a total derivative, so Qd reduces to 



a pure boundary term |28], which becomes negligible if the averaging is performed over 
a volume of the order of the Hubble volume, as should be the case if one is interested 
in the back-reaction of inhomogeneities on the background expansion rate. It was shown 



in Ref. [25| that in order for the back-reaction of cosmological inhomogeneities to drive 
acceleration, the kinematical back-reaction should be Qx> > 4:TtG{q)t>- 

There is another important feature of Eq. ( 5.11| ) which is worth pointing out, namely 



that a factor multiplies the mean spatial curvature term, implying that the Newtonian 
kinematical back-reaction Qxi couples with the PN mean spatial curvature {^^^R)^'^ , the 
PN kinematical back-reaction couples with the post-post-Newtonian mean spatial 

curvature {^^^R)^^^ and so on. 

By performing the conformal rescaling of the metric hap, the velocity-gradient tensor 
splits into a FRW and a peculiar term (recall ( |2.10D ) 

= i (^^ + ^5^) , (5.12) 

where a is the scale factor of the underlying FRW background and we use here the conformal 
time coordinate. In terms of the peculiar quantities, the kinematical backreaction becomes 

= 3^ ((^')^ - - ■ (5.13) 



Explicitly, for our one-dimensional Newtonian background we have a'^ =1} /3, thus Eq. (|5l^ ) 



becomes 

Si^ = -^Wl), (5.14) 

where a is our Einstein-de Sitter scale-factor, ?? = —Tdfip/{3{1 — T'^dfip/G)) and (. . .)t> is 
our Newtonian average 

(...)c = ^ / ...(1-t29iV/6) A, (5.15) 
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where Vv = J.j~, {l — d^q. Recalling that 



a 

where are proper Eulerian coordinates 
peculiar velocity in the direction ri, we can re-write equation (5.14) in the form 



(5.16) 

axa, vi = —Tdnp/3 is the Newtonian 



Qv 



(5.17) 



where we introduced the Eulerian average (. . .)^ = (1/ d^r) Jjy ■ ■ ■ d^r. Note that the last 
equality follows from the fact that {dr^vi)^ = 0, by construction (see Ref. @|). It follows 
that any analysis of back-reaction based on the Newtonian approximation is irrelevant. 
In addition, first and second order calculations based on the fully relativistic perturbation 
theory are not sufficient for evaluating back-reaction. This means that fully non-linear sub- 
Hubble perturbations must be taken into account. For these reasons, we use our non-linear 
PN metric to estimate the back-reaction term ( |5.13| ) in the PN approximation. 

In order to calculate the PN kinematical back-reaction Q^^, we consider the kinemat- 
ical quantities related to the PN peculiar velocity-gradient tensor 



1 



(5.18) 



In full generality, without any assumption about the Newtonian background, the PN ex- 
pansion of the average integrals in ( ^.13 ) leads to 



(5.19) 



1 



1 



1- 



1 



where with (. . we indicate the Newtonian average (PN corrections in the averaging 
procedure would yield higher-order terms). 



^2 



Recalling that in our case a = {) /3, equation (5.19) reduces to 



2 ,^ 



2 — q/ 



where we used the commutation rule for w 

{wYt, - {w')v = {wl)v - {w)vmv 



(5.20) 



(5.21) 



and the last equality in Eq. (5.20) follows from the previous Newtonian calculation, having 
neglected the last term, which involves {i!^)v- Finally, using our Newtonian average ( 5.15| ) 

,2' 
2 ) 

10ar2 



and our solution for Wo , we find 



8lVv Jv 



diidiipfd^q 



(5.22) 
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We can alternatively express this result in terms of Newtonian peculiar velocity vi. We 
find 

Qr = ^R,^?>i, (5.23) 

which indeed indicates a negligible PN correction. 
Finally, using the PN spatial curvature 

the PN expansion of Eq. ( ^ ) leads to 

^iS^PN^^f Q2^^3^_ (5.25) 

It is then clear that, in the case of planar dynamics, both the kinematical back-reaction 
and the mean spatial curvature reduce to purely boundary terms ^ and, as a consequence, 
they cannot lead to acceleration. Also interesting is the fact that {^^^R)^^ oc a^^, which 
provides a self-consistent solution of the integrability condition Eq. ( ^.llD , for vanishing 
Qv- 

5.1 Post-Newtonian corrections to the average scale-factor and expansion rate 

An alternative approach to back-reaction is that of computing directly the PN contribu- 
tion to the average scale-factor and expansion rate. Indeed, using the definition a© = 
(Vv/Vvo)^^^ , we can calculate the PN expansion of the average scale-factor, written as 

a-v = av + \aT>^ ■ (5.26) 



For the Newtonian term, we have for the volume 

^2 

.3 



Vv = a-' (^1 - -dtvj d'q , (5.27) 
thus, neglecting the boundary term, we obtain 

= — , (5.28) 



^0 



which is just the Einstein de Sitter scale factor a(r) 
For the PN term we find 



a^^ = ia(Hi,-(z/;o)7?o) > (5-29) 



where (. . is the Newtonian average and our solution gives 

-) + 1) (9iV 
6 — T^dfip 



63 ' "^IV^ y^^i^' ^ 3 

W = — 



fr^^iV (di^f + f (4 (ajvL - 1) + 1) (^K^)' + (4 (ajvL - 1) + 2) ifdfip 



(5.30) 



4 



Note that these quantities can be seen as boundary terms both in Lagrangian and Eulerian coordinates, 



since 



at the Newtonian level - [dF[q\)/dq )dq = (dF{qi{xi,T)/dx)dx, for any function -F(gi) 
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Explicitly, we have 



,4 



+ (4 (a^vL - 1) + 2) vj^iV] I , 
which can be written as 

up to negligible boundary terms. The average scale-factor then becomes 

av = a{l--^{{l+I)vl)v^, (5.33) 
which indicates a negligible PN correction. 

The PN contribution to the average expansion rate (G)x) = {3/a)a'jy/ax> is obtained by the 
expansion of its very definition, leading to 

{@)v = 3H + ^{wyj,, (5.34) 

where H{t) = 2/(ot) is the Hubble expansion rate of our Einstein-de Sitter background. 
Neglecting boundary terms, it is straightforward to obtain ^ 

(G)^ = 3h(^1- -^{{l + 6)vi)v^ . (5.35) 

This result again implies that the PN correction is fully negligible for plane parallel per- 
turbations. 

6. Comparison with the Szekeres solution 

In this paper we have considered the evolution of an irrotational and collisionless fluid 
in General Relativity in the synchronous and comoving gauge. Following the fluid-flow 
approach, it is possible to alternatively describe our system in terms of the fluid 

properties of irrotational dust, i.e. mass density, volume-expansion and shear tensor, and 
the electric and magnetic parts of the Weyl tensor. In addition, in the special case of plane- 
parallel dynamics considered here, the magnetic Weyl tensor vanishes identically, thus 
leading to the so-called silent universe case [32, 33, 34 1, described by the set of equations 



Q = -Qq (6.1) 

e = -QT? - \q'^ - AnGg (6.2) 

s = i;2 - ^es - E (6.3) 

E = -3ET. -@E- A-kGqT. , (6.4) 



^If we expand Eq. (5.35) up to second-order in perturbation theory, we find quantitative agreement with 
Eq. (41) inRef. We checked that the difTerent numerical factor is simply due to the different definition 
of {S9)t> used in the two calculations. 
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where is the trace of velocity-gradient tensor, S and E are the eigenvalues of the shear 
tensor and of the electric Weyl tensor 



= H% - e') + ee^ - e^e} + c'( - ^'^r ) (6.5) 



in the directions q2 and q^. In this framework Croudace et al. |35| obtain what they refere 
to as relativistic Zel'dovich solution, a sub-case of the exact solutions by Szekeres, p6|| . 
This solution is made more appealing by the recent growing interest on the Szekeres metric 
in the general framework of studying inhomogeneous cosmologies as possible alternatives 



to FRW (see, e.g. Refs. [^] and |35] and Refs. therein). Ref. [35] considers the relativistic 
evolution equations of silent universes. In the special case of local planar symmetry, i.e. 
^2 = A3 = 0, where Aj are the eigenvalues of the matrix dad pip, the velocity-gradient tensor 
reads || 

1 



gAi 
1— aAi 



e^ = -| 1 I, (6.6) 



a 



1, 



which corresponds to the Newtonian Zel'dovich solution. Croudace et el. then compute 
the associated metric from 



1 



i/3<7 (6.7) 



2 

via the time-time component of Einstein's equations, i.e. the energy constraint, that closes 
the relativistic fluid- flow equations, completely fixing the spatial dependence of the metric. 
They show that the resulting metric coincides with the Szekeres form 

h22 = /i33 = 1 (6.8) 

hii = {d{qa) - a{t)c{qi)f (6.9) 

with 

d{qa) = din((?l) - \c{qi){{q2f + (qsf) - a{t)c{q{) . (6.10) 

This solution leaves the FRW expansion unperturbed in the directions q2 and 53. In 
fact, following Ref. |Q, Croudace et al. discard the sub- leading PN trace part of the 
initial conditions ( |4.1C| ) , thereby considering perturbations in the direction qi only. On the 
contrary, we kept the one-dimensional initial seed ip{qi) in all directions, thereby allowing 
for perturbations in the component /i22 and /i33 of the metric ah initio. As a consequence, 
our solution for /122 and /i33 changes with time, showing a non-linear dependence on the 
gravitational potential ip. Hence our (approximate) solution cannot be recast into the 
Szekeres form. 



7. Conclusions 



In this paper we have obtained a solution of Einstein's field equations describing the non- 
linear cosmological dynamics of irrotational dust in Lagrangian coordinates for the specific 
case of plane-parallel perturbations, which is exact up to first post-Newtonian order. Our 



solution (Eq. (|4.37] )), which represents the post-Newtonian extension of the well-known 
Zel'dovich solution for plane-parallel cosmological dynamics in Newtonian gravity, can be 
useful in many respects. Here we have seen just one possible application, in the frame of 
averaging and back-reaction of cosmological inhomogeneities. In a forthcoming paper we 
will analyze the properties of photon geodesies in our metric, to study, for instance, the 
luminosity distance-redshift relation in a strongly inhomogeneous Universe such as the one 
described by our metric. 

Our results also allow us to study the final stages of plane-parallel collapse, which ob- 
viously leads to a shell-crossing, pancake-like singularity. Caustic formation is considered 
the main limitation of the Lagrangian description, both in the Newtonian approximation 
and in General Relativity. As well-known, caustics arise because several fluid elements 
coming from different positions may converge to the same Eulerian position, thus forming 
infinite density regions. Such a pathological behavior occurs in our case when / = — 1, 
i.e. when r^c?^(^/6 = 1: at this time the determinants of both the Newtonian (4.6) and 
PN metric ( [4.37| ) go to zero, while the density contrast at both the Newtonian, Eq. ([4.^), 
and PN level, ( [4.26| ), becomes infinite, as does the PN spatial curvature in Eq. ( [5.24| ). The 
appearance of shell-crossing singularities can be understood as indicating the breakdown 
of the dust approximation, rather than the occurrence of a true physical singularity of the 
gravitational collapse. The formation of caustics appears as an artifact of the extrapolation 
of the pressure-less fluid approximation beyond the point at which pressure has become 
important. An important result of our analysis is that the divergence of the PN spatial cur- 
vature at caustic formation is completely eliminated by the spatial smoothing procedure. 
Indeed TZ^^ diverges like {1 + 5), which is exactly compensated by the square root of the 
spatial metric determinant oc 1/(1 + 5). This very fact confirms that the instability found 
in Ref. using a gradient expansion technique, and in Ref. using a different ap- 
proximation scheme, cannot be interpreted as a consequence of a shell-crossing singularity, 
but really arises from the back-reaction of sub-Hubble modes. 

Our results here bring both good and bad news for the back-reaction of cosmic inho- 
mogeneities to represent a potentially viable alternative to dark energy at a fundamental 
level. Let us start with the bad news: as we have shown the kinematical back-reaction 
scalar remains a negligible boundary term also at the post-Newtonian level, so that no rel- 
evant back-reaction effect is implied by our PN solution. We should stress, however, that 
such a result is clearly a consequence of our very specific (but analytically solvable) model, 
which relies on one-dimensional dynamics. It is a very reasonable hypothesis that such a 
result will be modified by a (more complex!) full 3D calculation; this will be the subject 
of a subsequent paper. The good news are: the very fact that the Lagrangian approach 
allows to obtain a quantitative estimate of back-reaction and the fact that the appear- 
ance of shell-crossing singularities does not lead to divergences in the average dynamics of 
inhomogeneous dust Universes. 

Also interesting, in our opinion, is the fact that the PN metric depends on the pri- 
mordial non-Gaussianity strength parameter cnl (which is related to the best known non- 
Gaussianity parameter /nl by /nl = (5/3)aNL! for large onl)- Although in our final results 
for , {^^^R)^^ and this quantity disappeared (in that it only enters in boundary 
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terms), one can safely predict that this parameter wih affect back-reaction quantities for 
general three-dimensional perturbations, hence suggesting the intriguing possibility of a 
potential quantitative link between primordial non-Gaussianity and back-reaction. 
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